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Abstract. This paper deals with a class of macroscopic models for cell mi- 
gration in a saturated medium for two-species mixtures. Those species tend to 
achieve some motion according to a desired velocity, and congestion forces them 
to adapt their velocity. This adaptation is modelled by a correction velocity 
which is chosen minimal in a least-square sense. We are especially interested 
in two situations: a single active species moves in a passive matrix (cell migra- 
tion) with a given desired velocity, and a closed-loop Keller-Segel type model, 
where the desired velocity is the gradient of a self-emitted chemoattractant. 

We propose a theoretical framework for the open-loop model (desired ve- 
locities are defined as gradients of given functions) based on a formulation in 
the form of a gradient fiow in the Wasserstein space. We propose a numerical 
strategy to discretize the model, and illustrate its behaviour in the case of a 
prescribed velocity, and for the saturated Keller-Segel model. 

1. Introduction, modelling aspects. We propose a model to handle interactions 
between organisms such as unicellular organisms e.g. bacteria or amoebia. The be- 
havior of an entity is based on a will to move with its desired velocity, regardless of 
others, but the fulfillment of individual wills is made impossible because of conges- 
tion. We consider here the case of two organisms (see Section 4 for a straightforward 
extension to several populations). We describe the densities by measures pi and p2. 
We assume global saturation of the domain, i.e. pi + P2 = 1- In order to take into 
account the individual tendencies together with congestion constraints, we consider 
that individuals have two velocities, namely a desired velocity that will be denoted 
by Ui {resp. U2) and a common correction velocity that will be denoted by w. 
Densities pi and p2 satisfy: 

atP^ + V-(p,(U^ + w)) = z = l, 2. (1) 

We consider the correction velocity which minimizes norm among all those ve- 
locity fields which ensure preservation of the saturation constraint pi -\- P2 = 1, 
which can be expressed in a dual, Darcy-like, form: 

r w + vp = 

\ V-w = -V • (piUi + . ^ ' 
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This basic model corresponds to a competition between two species which tend 
to achieve a prescribed motion, and reahze a sort of compromise. We shall be 
especially interested in the following situations: 

(i) Competition between species which tend to minimize a given function: is 
defined as — Vl^i, with piUi + P2U2 not feasible (i.e. leading to violation of 
the saturation constraint) in general. This setting, which extends the model 
proposed in [11] to a two-population situation with distinct tendencies, shall 
be the core of the theoretical analysis proposed in section 2. 

(ii) Desired velocity for species 2 is zero, and Ui is prescribed as the gradient of 
a given function S (concentration of a chemoattractant). Note that this case, 
which corresponds to the migration of cells in a passive biological matrix, as 
encountered in the developpment of atherosclerosis, is a particular case of {i). 

(iii) Keller-Segel model with congestion (closed- loop version of the basic model): 
U2 is again 0, and Ui is defined as V^*, where 5* is a chemotactic agent created 
by the species 1 itself. Assuming S diffuses in the mixture domain, we shall 
consider 

dtS-AS = pi, 

or, assuming diffusion is instantaneous, — A^* = pi. 

The question of boundary conditions rises delicate issues. Setting equation in 
the whole space amounts to deal with infinite quantities of 1 and/or 2 (the sum of 
densities is equal to 1), which rules out standard tools for theoretical analysis and 
numerical simulation. One may consider that both populations occupy a bounded, 
moving domain Q{t), on the boundary of which pressure is zero. It is then natural 
to consider that the boundary of Q moves with a normal velocity which identifies 
to the normal velocity of the mixture (piUi + P2U2 + w) • n. This approach is 
surely relevant in some situations, yet we shall not consider it in this paper. Let 
us give an idea of one of the difficulties which are likely to occur: consider the 
monodimensional situation with initial condition and desired velocities: 

Pi = 1(0,1) , P2 = l(-i,o) , Ui = 1 , U2 = 0. 

Translation of those intervals at constant speed 1/2 is obviously a solution to our 
problem (with piecewise affine pressure, at both ends, 1/2 at the interface). Now 
consider the scenario where 1 runs away with velocity 1, and 2 stays where it is 
(both species are separated, and the domain Q{t) is no longer convex). One can 
even imagine that half of 1 (initially located in 1 (1/2,1)) takes off with velocity 1, 
while the other half of 1 trails 2 at speed 1/3. In this spirit, one can create an 
infinite number of scenarios which all seem to be, in a reasonable sense, solutions to 
our problem. Note that those problems are likely to occur in the case of expansion 
fields, as the one considered previously, because of the loss of monotonicity of the 
underlying evolution equation. 

Another approach consists in considering a fixed domain delimited by rigid 
walls. This assumption is reasonable if one considers migration phenomena in living 
organisms. The boundary condition for Darcy equations is of Neumann type: no- 
fiux condition at walls leads to (piUi + P2U2 + w) • n = 0, i.e. 

9P / XT XT N 

— = (piUi+p2U2)-n. 

This condition ensures no-fiux for the mixture, but not for both species indepen- 
dently; it may allow some out- or in-fiux of either 1 or 2, which calls for further 
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prescriptions on the boundary. The situations we shall consider in next sections 
suggest that, at least in the case U2 = 0, the model tends to create pure zones (ei- 
ther 1 only, or no 1 at all) in the neighbourhood of boundaries. As a consequence, 
global no flux conditions turn into no flux conditions for both species. Yet, in 2 or 3 
dimensions and for non zero desired velocities, one may have to consider more gen- 
eral situations. To overcome those difficulties, we shall consider in the theoretical 
part (Section 2) the fully periodic setting. 

In a different context, such saturated two-component mixtures have been studied 
in [13], with a motion driven by chemical potentials. In the biological context, a 
Keller-Segel model with logistic sensitivity (KSLS) has been proposed recently [6, 4]. 
This KSLS model reads 

dtp + V-{p{l-p)\J) = 0, (3) 

where the "desired" velocity U is of the chemotactic type: 

U = V5, dtS-AS = p. 

First of all, in one dimension (say = (0, 1)) with no-flux boundary conditions 
at X = 1, our model reads 

-dxxP = -dxipU) 

with (no flux condition) pU — dxP = at 1, so that dxP = pU in the whole interval, 
and therefore it identifies exactly with KSLS. Note that in case the velocity is given 
and constant, we recover an inviscid Burgers equation. 

Both models are different for d > 2. Yet, they present some common structure 
which can be described as follows: consider (situation (Hi) above) the case of zero 
velocity for species 2, and velocity for 1 given as the gradient of a density S of 
chemoattractant created by 1 itself. We consider the biperiodic setting, and de- 
fine as the operator which maps a function g with zero mean value over the 
biperiodic domain onto the zero mean value solution to 

Au = g. 

The pressure in our model can then be written p = V • (/>U), so that transport 
equation for species 1 becomes 

StP + V- {p (U- VA-iV-(pU))) =0. 

which is formally equivalent to KSLS model (3) with identity operator replaced by 
Helmoltz projection VA~^V- (projection onto the space of irrotational fields). Non- 
locality of the latter operator differentiate both models. Yet both PDE's systems 
are quite similar from the spectral point of view. 

As for modelling aspects, they rely on different assumptions: KSLS model can 
be seen as standard KS model with a reduced velocity U = (1 — p)V5', which can 
express for example the fact that when entities reach some critical density (set to 
1 here), they disturb each other and are no longer able to sense the underlying 
gradient. In our approach, we consider that entities continue to exert some action 
which would lead them in the right direction if they were alone, but they are (me- 
chanically) prevented from fulfilling their purpose, because of the presence of other 
entities (possibly with different "strategies"). The model is in some way less con- 
strained as motion of saturated zones is possible, but also more constrained because 
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the other species (which replaces empty space in KSLS) has to be swept away, at 
some price, by species 1. 

2. Theoretical framework. We consider in this section the case of 2 species in 
the flat torus ft = x S^, with desired velocities prescribed as follows 

Ui = -VI^i, U2 = -VI)2, 

where Di and D2 are given smooth functions over Q. The system corresponding to 
this situation writes 



f StPi + V.(pi(Ui+w)) = 

StP2+V-(p2(U2+w)) = 



-Ap = -V-(piUi+p2U2). 

Although concentrations for both species automatically admit a density which 
is in (thanks to the congestion constraint), we shall keep the general setting 
of measures, which is usually adopted in optimal transport, and which allows for 
generalisations (e.g. relaxing the congestion constraint in some zones). Besides, 
we disregard normalization to obtain probability measures: we shall consider that 
pi 0p2 belongs toV{^xQ)^ although the total mass is not 1. Similarly, we consider 
that Pi e V(Vt). 

In what follows, regularity of functions defined over Vt accounts for periodic 
conditions. In particular H^(Q) is defined as the set of biperiodic functions with 
regularity (closure for the norm of biperiodic regular functions). 

Definition 2.1 (Weak solutions). We say that (pi,p2) is a weak solution of (4) 
with initial condition (/>i,/)2) if Pi + P2 = 1 for a.e. t, and if there exists p G H^(Q) 
such that for ah Lpi,Lp2 G C^([0, T[xl]), for a.e. t G (0,T), and for ah q G H^{^), 
we have 

/ / {dt^i + V^i-{\5i-Wp))dpi + [ M^,.)dp\ = 
Jo Jii Jn 

f I {dm + V^2 ■ (U2 - Vp)) dp2 + I ^2(0, .)dpl = 0, 

1 Vp'Vq = I Vq-Uidpi ^ / Vq'\J2dp2. 
Jn Jn Jq 

Let p = pi^p2 for a.e. t, U = (Ui, U2). We can rewrite the previous formulation 
as follows: for ah (p G C^{[0,T[x{n x Q)) 

r I {dt^ + (U - (Vp, Vp), V^)) dp ^ I (^(0, .) dp"" = 0. 
Jo JQxQ JnxQ 

and for a.e. t G (0,T), ah q G H^{n), 

I Vp Vq = I ((Vg,Vg),U)dp. 

JQ. JQxQ 

Let us state the main result of this section: 
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Theorem 2.2 (Existence of a solution). If Di and D2 are Lips chit z functions, and 
{pi-> P2) positive measures that satisfy p\-\- p2 = l a.e., then equation (4) admits at 
least a weak solution, according to definition 2.1. 

The proof of this theorem rehes on the notion of gradient flow in the Wasserstein 
space (see e.g. [1]). Theoretical analysis of a similar problem in the context of crowd 
motions was proposed in [11]. The proof proposed therein is quite general, and the 
same approach could be carried out in the present context. Yet, to avoid some 
technicalities and to give a clearer view of the underlying gradient flow structure, 
we propose here an alternative approach which is directly based on the main exis- 
tence and characterization theorems in [1]. The rest of this section describes the 
main outlines of this approach, which is based on theories of optimal transport and 
gradient flows, which we will use to prove theorem 2.2. For more details, see the 
books of Villani [17], [18], and Ambrosio et al. [1], [2]. See also the recent appli- 
cation of this framework [7] to recover entropic solutions to the Burgers' equation 
(which identifies in some way to the proposed model for (i = 1, as pointed out in 
the introduction). 

We define a discrete scheme for (4): 

Definition 2.3 (JKO scheme). Let r > be given, and 

an initial density. We define , . . . , p"^ (we drop the dependence upon r to alleviate 
notations) recursively according to 

p^e argmin f J(p) + /^(p) + -1 VP^|(p, p-^)^ (5) 



where J is given by 



J{p) = / {Di{xi) ^ D2{x2))dp{xi,X2) 

and K is the set of admissible densities 

K = {p e X ft) : p = pi (S) P2j Ml + = 1 a.e.}. 
The notation Ik stands for the indicatrix function of i.e. 



if peK 

foo if p ^ K. 



This is a well-known scheme in gradient flow theory (see [5], [8], [3], [1], [2]), which 
has been widely used to prove existence theorems. Taking a minimizing sequence 
of the right-hand side of (5), we can verify that every of these minimizing problems 
admit at least a solution. 

Let us underline that here the JKO scheme is only a tool to prove the existence 
of a solution to our problem, and has not been used for numerical purposes. As 
a matter of fact, the numerical resolution at each time step of the minimisation 
problem leads to many difficulties, and we have chosen a totally different approach 
for the numerical tests in section 3. 

Proposition 1 (Distance between two product measures). If p = pi <^ p2^i^ = 
z^i ^1^2, then 

Wi{p,v) = Wi{pi,yi)^Wi{p2^V2)- 
Moreover, if yi (resp. V2) is the optimal transport between pi and vi (resp. p2 and 
V2), then r = (ri,r2) is the optimal transport between p and v. 
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Proof. As (ri, r2) transports fi to the previous definition of W2 gives W2^(/i, 1^) < 
1^2 (Mi 5 ^1) + ^2 (M2, ^2)- To obtain the converse inequality, we simply use the dual 
definition of W2 with Kantorovich potentials (see [17]). □ 

Thanks to the previous proposition, the JKO scheme can be rewritten as follows 

(p?,pj)e argmin ( [ D,dp, + [ D2dp2 + ^WHp^, p^') + ^Wiip2, pr')) 
Pi+P2=i Vio Jq It J 

(6) 

Proposition 2 (Convergence of the JKO scheme). Let pr be the piecewise constant 
interpolation (in time) of the sequence {p^)n defined by the JKO scheme{5) for the 
time step r. Under the assumptions of theorem 2.2, there exists a subsequence of 
{Pt)t which converges weakly to a weak solution of the transport equation 

dtP^V - (pu) = 

where u verifies 

ue-d{J^lK){p). 



Proof. We have to verify that the functionnal ^ := J-\-Ik satisfies the assumptions 
of the theory developed by Ambrosio et al. in [1]. ^ is clearly proper and lower 
semicontinuous, and its sublevels [^{p) < c] are compact. Moreover, we can prove 
that $ is regular according to definition 10.1.4 in [1], i.e. for all sequence = 
(Pi P2) ^ ^5 for all G d^{p^) such that (p^) narrowly converges towards 
p in V{ft X 17), and 

sup ||u^||x,2(^^) < +00, ^ u weakly 

n 

then u e d^{p). All assumptions of Prop. 2.2.3, Th. 2.3.1, and Th. 11.1.3 of [1] are 

satisfied, therefore there exists a subsequence of (pr)r which converges to a weak 
solution of 

atp + v-(pu) = 0, 

with 

u G —d^{p) for a.e. t. 

□ 

Proposition 3. The limit of {Pt)t proposition 2 is actually a weak solution of 
equation (4). 

Proof. It is quite easy to verify that when r converges to 0, the properties p(t, .) = 
Pi{t,') ^ P2{t,') with pi(t,x) ^ p2(t,x) = 1 for a.e. x, and \i{x,y) = (ui(x), U2(^)) 
still hold true. Therefore, the continuity equation rewrites 

/ S,pi+V-(piUi) = , . 

1 atP2+V-(p2U2) = ' 

and the condition p^ + p2 = 1 gives the following equation on ui and U2 

V • (piui +P2U2) = 0. 

We now use the fact that — u is a strong subdifferential, i.e. that for all transport 
map r, we have 

$(p)- / (u,r-id)p < <l>(r#p) + o(||r-id||i,2(^)). 
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Let us underline that if r^p the previous inequahty does not give any infor- 
mation, as <l>(r^p) = +00. 

We define the set of admissible velocities as follows 

Cp := {v = (vi,V2) G L'^ipi) X L'^{p2) : V • (piVi + P2V2) = 0}. 

We just proved that u e Cp. Let v G Cp, £ 7^ 0, and = id + ev. In most 
cases, we have r^p K. However, it is possible to find a transport such that 
(t^ o Y^)^p G and W2{{t^ o T^)^p^T^p) = o{e). The subdifferential inequality 
applied to o gives 

/ {Di + D2) dp- (u, o - id) dp 

< I (A + Z?2)(t^or^)#^ + o(||t^or^-id|U.(,)). 
Using estimates between (t^ ov^)^p and r^p, we get 

/ {Di^D2)dp- [ {u,r' -id) dp < [ {Di^D2)v%p^o{e), 
jQxn JQxn Jnxn 

which implies 

-\/D -u,ev)dp < o{e). 



In: 



/nxQ 

Letting e go to from positive and negative values, we obtain 



/ {-\/D -u,v)dp = 0, 



which means that —VD — u G C^. As u G Cp, u is indeed the projection of 
-VD = (Ui,U2) onto Cp, i.e. 



u = argmin / |U — v| dp. 

This minimizing problem can be written as a saddle-point problem for which we 
can prove existence and uniqueness of a solution (u,p) which satisfies 

u + (Vp,Vp) = U, 

i.e. 

r ui = Ui - vp 

[ U2 = U2 - Vp 

Moreover, since V • (piUi + P2U2) = 0, the equation on p is given by 
-Ap = -V • (piVp + p2Vp) = -V • (piUi + P2U2). 
If we replace the expression of u in equations (7), we finally get 

StPi+V.(pi(Ui-Vp)) = 

StP2+V-(p2(U2- Vp)) = , 

-Ap = -V- (piUi +P2U2) 
i.e. p is a solution of (4). □ 
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Remark 1. (Uniqueness) We focused here on the proof of existence of a solution. 
Under reasonable assumptions on Di and 1^2, it is to be expected that the JKO 
process leads to a unique solution (see [11] for remarks regarding uniqueness in a 
similar setting). Yet, as we pointed out in the introduction, the system is under 
some conditions equivalent to the standard inviscid Burgers' equation, which rules 
out uniqueness in general. It is natural to wonder whether the JKO scheme selects 
a particular solution, namely the entropic one. This delicate question is still widely 
open, but a recent work on Burgers equation ([7]) suggests a positive answer. 

3. Numerical methods and results. We consider here the transport model with 
congestion, with one active species only (Ui = U, U2 = 0) expressed in terms of 
the active species only: 

dtp + V ■{p{\J + w)) = 0, (8) 
w = -Vp, (9) 
-Ap = -V-(pU). (10) 

3.1. Discretization and maximum principle. First, a time discretization of the 
system (8)-(10) is given. We set t" = nSt, and 

p'^ix)^ p{e,x), w"(a;) ~ w(t",a;), U"(a;) ~ U(r, x) , p"(a;) ~ a;) . 

Equation (8) is approaciied by using a backwards Euler finite difference method. 



Semi-discretized version of System (8)-(10) hence reads: 

pu+i ^ ^"-,5tV-(p"(U" + w")), (11) 

w" = -Vp", (12) 

-Ap" = -V-(p"U"). (13) 



For the sake of simplicity, we only give the ID spatial discretization of (11)- (13), 
the extension to 2D or 3D on cartesian grids being straightforward. 
We introduce: 

Since the equations of the model are written in a conservative form, the natural 
framework to be used for the spatial discretization of (11)- (13) is the finite volume 
framework. We hence introduce the control volume defined by: 

We denote by (p^)i and {p2)i the piecewise constant approximations of densities 
and pressures at time (i.e. stands for the value of p at cell center x^, at time 
t^). As for fiux variables, we define approximate values at interfaces: (resp. 
U^) stands for w (resp. U) at interface x^_i^ 

By integrating (11) over d and by using the divergence theorem, we obtain 

/ p^+Mx= / p^dx-^t [p^(U^ + w^)]I^^* . (14) 
We approximate the integral terms in the following way: 

JCi JCi 
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An upwind discrete flux is then introduced to approximate the remaining flux term 
in (14): 

p-{x,_.) (u-(x,_.) + w-(x,_i)) ^ A^^ (Ur,pr-i,pr) (wr,P?-i,P?) 

(15) 

where the numerical flux A^^ is defined by 



A^P{u, p-, p+) = 

We finally obtain ; 




if > , 
ifu<0. 



- Yx (w?+i,P?,P?+i) - (w?,p?_i,p?)) . (16) 

It is important to notice that fluxes pU and pw were treated separately in (15). 
As we will see, this plays an essential role in preserving the maximum principle on 
p for the numerical solution. This advection scheme is stable under the Courant- 
Friedrichs-Lewy condition : 

8t<\- , ^ — . (17) 

2 |U«|oo + |w«|oo ^ ' 

Eq. (13) is discretized in space with 

p-Vi - 2g + ^ ^ (ur+i, pr, pr+i) - (ur, ^r-i, ^D) • m 

Finally, by using an Euler finite difference scheme in (12), correction velocity w is 
approximated with 



Proposition 4. The numerical scheme (16)-(19) satisfies the following maximum 
principle: 



0<p?<l, Vz = l,...,iV, 
then < < 1, Vi = 1,...,7V^ Vn. 

Proof. For the sake of simplicity, only periodic boundary conditions are considered 
in the following proof: 

Pi = Pn. > 
Pi=Pn.- 

Note that, in practical, the proposition remains true if the no-flux boundary con- 
dition (mentioned in Section 1) is taken over p. In this case, appropriate boundary 
conditions have to be taken upon p in the upwind scheme: 

A-P (iii,Po,Pi) = 0, 

^""^ (^iVx + l,PA/-.,PA/-. + l) =0- 

First the positivity of the numerical scheme is given by the well known positivity of 
the upwind scheme under the C.F.L. condition (17). For the upper bound p^ < 1, 
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the proof rehes on the positivity of pf := 1 — that is given by a similar numerical 
scheme. Indeed, equation (16) implies 

+ ^ {A-p 1 - M?, 1 - A^r+i) - (w?, 1 - ^ti, 1 - P?)) ■ 

Since A'^p{u^ 1 — — p~) = ii — A^p{u^ p'^,p~) we have 

- £ (^"^ (wiVi,M?,M?+i) - (wr,M?-i,Air)) • 

By replacing the values of with the discrete gradient of p as shown in (19), one 
gets 

p7^' = + S (^"' (ur+i,p?,pr+i) - A-p (ur,p?_i,pr)) 

St -2pf 



(5a; \ i5a; 

- ^ (^"' (wIVi,Ai?,Air+i) - A-P (w?, • (20) 

We recognize in the first terms of the right hand side of (20) the exact discrete 
expression of (18). Hence the numerical scheme satisfied by p is 

p7^' = - ^ (^"' - «,P?-i,P?)) , 

which is an upwind discretization of the advection problem: dtp -\- V ■ (pw) = 0. 
Thanks to the positivity of the upwind scheme under the C.F.L. condition (17) we 
deduce the discrete maximum principle for p. □ 

3.2. Numerical results. In this section we present several numerical simulations 
performed for the migration model (8)- (10) with the numerical scheme (16)- (19) 
introduced above. First, in order to validate the numerical method, two ID test 
cases for which the exact solution is known are presented. Then 2D simulations 
are shown for two types of situations: first the case where the desired velocity U is 
known, and then the case where U is given as a function of p. 

3.2.1. ID simulations. In this section the domain is the interval (0, 1), the desired 
velocity is U = 1, and we consider the following boundary conditions for p: 

p^(x = 0)=0, 5,p^(x = l)=p^(l)U^(l). (21) 

Two different initial conditions are considered: 

= ^l[o.i,o.9] +l[o.9,i] (see Fig. la) , (22) 

and 

p' = l[o.3,o.5] (see Fig. 2a). (23) 

As we can see on Figures Ib-c and 2b-c, errors between the computed and the 
exact solution happen mostly near the points of discontinuity of the exact solution, 
except when we reach steady state. Discontinuities of steady states are indeed 
captured with striking accuracy (see Figs Id and 2d), whereas upwinding might be 
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Figure 1. Evolution of the concentration p for the initial data 
defined in (22), with the "wall" boundary condition (21). Exact 
solution (plain) and discrete solution (dotted hue). Simulations 
performed for Nx = 200. 



expected to diffuse fronts: it is due to the behavior of the model we consider, which 
tends to sweep p toward the discontinuity, inducing permanent correction of the 
numerical diffusion. 

3.2.2. 2D simulations for a constant desired velocity U. In this section we perform 
some numerical simulations for the migration model (8)- (10), in the case where the 
desired velocity U is given and constant in time. In all the cases shown here, the 
domain Q is bounded and surrounded by walls, so that p verifies Neumann boundary 
conditions: 

| = ,U.n. (24) 

The first situation studied here is a direct 2D extension of the ID test case described 
previously, see Fig. 2. The domain is the unit square 1^ = (0 , 1) x (0 , 1). We also 
prescribe the following desired velocity: U = 1. Figure 3 shows the evolution of p 
from the initial condition (Fig. 3a) to the steady state (Fig. 3f) given by a numerical 
simulation performed in the situation described above. The same way as in the ID 
situation (Fig. 2b), the evolution of p begins with the spreading of a mixing zone 
along the x axis (Figs. 3b-c). Along with the evolution of the mixing zone, the 
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Figure 2. Evolution of the concentration p for the initial data 
defined in (23), with the "wall" boundary condition (21). Exact 
solution (plain) and discrete solution (dotted hue). Simulations 
performed for Nx = 200. 

solution also tends to spread along the y axis. Hence, a real difference is shown here 
between our model and Burgers-type models. When the mixing zone reaches the 
wall, a congested zone instantly appears and develops backwards (Fig. 3d-e). The 
steady state is reached when all the right side of the domain is saturated with the 
active species (Fig. 3f). 

The second situation studied in this section involves a more complex geometry 
for the domain, which is defined by 

1] = (0 , 1) X (0 , 1)\ ([0.3 , 0.7] X [0 , 0.45] U [0.3 , 0.7] x [0.55 , 1]) , (25) 

which corresponds to two rectangular rooms joined by a thin corridor (see Fig. 4). 
The desired velocity is set to U = —VI), where D{x^y) represents the geodesic 
distance from (x, y) to the right wall. In practice D is computed thanks to the 
toolbox provided by [15], which is based on the Fast-Marching method (see [10]). 
Figure 5 represents the simulated evolution of the concentration p with the setting 
described above. The initial condition is shown in figure 5a. As we can see a 
congested zone rapidly forms near the entrance of the corridor (figure 5b), which 
size decreases in time as the matter runs through it (figures 5c-e) . We notice that a 
quasi- homogeneous flow regime tends to be established in the corridor (Fig. 5d-e), 
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Figure 3. Evolution of the concentration p in the case of a con- 
stant desired velocity U = 1, with the "wall" boundary condi- 
tion (24). Mesh resolution: A^^ x A^^ = 300 x 300. 



where a constant concentration p ~ 0.5 is observed. It illustrates conservation of 
species 2 (in white in the figures): as the domain is bounded by walls, flow of 1 
to the right has to be balanced by an opposite flow of 2 to the left. As we define 
the correction velocity in the least-square sense, actual speeds for 1 and 2 are close, 
so that the mixture is necessarily balanced for global mass conservation reasons. 
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Figure 4. Geodesic distance D{x^y) to the right wah taking into 
account the obstacle computed with a fast marching algorithm. 



After some time, as observed in the previous test case, all the active species is 
concentrated on the right side of the domain (Fig. 5f). 

3.2.3. 2D simulations for an evolving desired velocity U(p). As mentionned in sec- 
tion 1, the desired velocity can be chosen to depend on the local concentration p, 
the same way as in the Keller-Segel model (see [9]): 

U = V5, (26) 

A5 = -p. (27) 

In this section we consider periodic boundary conditions over (p,p, 6'): 

p{t, 0, y) = p{t, 1, y) , p{t, X, 0) = p{t, X, 1) , (28) 

pit, 0, y) = p{t, 1, y) , p{t, X, 0) = p{t, X, 1) , (29) 

S{t, 0, y) = S{t, 1, y) , S{t, x, 0) = S{t, x, 1) . (30) 

The numerical treatment we give to equations (26)- (27) is similar to the numerical 
treatment of equations (9)- (10) described previsouly: 

qn en 



5x 

on _ qn 

on 9 on _|_ on on 9 on 1 on 

fe^ + " ^^^^ 

For all the following numerical simulations the initial data is set randomly according 
to Bernoulli's law 

=0) =1-9, (34) 

= !)=</, (35) 

where q G (0, 1) conditions the initial mass of active cells. The domain is the 
unit square, the boundary conditions set here are periodic boundary conditions 
for p, p and S. Figures 6 and 7 represent the evolution of p given by numerical 
simulations upon our model with the desired velocity defined by (26)- (27), and for 
initial conditions where g = 0.1 and g = 0.5 (respectively represented in Figures 6a 
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Figure 5. Evolution of the concentration p in the case of a con- 
stant desired velocity U = VI^ (see Fig. 4), with the "wall" bound- 
ary condition (24). Mesh resolution: N^^ x Ny = 300 x 300. 



and 7a) . In each case we see aggregation of bigger and bigger structures as time 
goes (Figs. 6b-e and 7b-e). Two different steady states are shown, both involving 
a steady, fully congested zone. In the case g = 0.1 the steady state congested zone 
is a disc, whereas in the case q = 0.5 a band has formed because of the periodic 
boundary conditions we chose to consider. 
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Figure 6. Evolution of the concentration p in the case of a desired 
velocity given by (26)- (27), with periodic boundary conditions (28)- 
(30). The domain is initially filled randomly with 10% moving 
species. Mesh resolution: Nx x Ny = 300 x 300. 



4. Conclusion, extensions. We proposed a model to describe the motion of mix- 
tures of cell populations in a saturated medium with a constraint on the local 
density. We provided an adapted theoretical framework, based on a reformulation 
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Figure 7. Evolution of the concentration p in the case of a desired 
velocity given by (26)- (27), with periodic boundary conditions (28)- 
(30). The domain is initially filled randomly with 50% moving 
species. Mesh resolution: Nx x Ny = 300 x 300. 



of the model as a gradient fiow in a product space of densities, and proposed a 
discretization strategy which enjoys reasonable stability and accuracy properties. 
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In terms of modelling, as mentioned in the introduction, any number N of species 
can be handled, as far as equations are considered (possible issues concerning bound- 
ary conditions are disregarded here). Saturation writes simply pi + • • • + Pat = 1, 
we have an advection equation for each species 

dtPi + V • (pi (U^ + w)) = , i = 1, . . . , TV, 

and the common correction velocity verifies 

V-w = -V - (jlp^^i^ ■ 

It could be of particular importance to include also proliferation phenomena. De- 
noting by Pi the local rate of creation (possibly depending explicity upon pi or other 
densities), conservation equations become 

dtPi + V • (p, (U, + w)) = ft , i = 1, . . . , AT, 
and the constraint on w accounts for mass creation: 

N / N \ 

i=l \i=l ) 

Note that, in this situation, global non conservation rules out the use of no-fiux 
conditions (or periodic setting). Considering a bounded domain ^ one can con- 
sider its boundary as a free outlet (interaction pressure is set at 0), so that a non 
conservative fiux through boundary can compensate the unbalance of mass in the 
domain. In case of mass creation (ft > 0), one can expect outflow through the 
boundary, so that no additional condition is needed for the advection equations. In 
case some individual velocities U^+w may point inward the domain, the system has 
to be complemented with appropriate conditions (e.g. prescribed value of ingoing 
densities). 

As for theoretical issues, non conservation rules out the standard framework of 
gradient flow in the Wasser stein space, which is dedicated to measures with constant 
total mass. Yet, as suggested in [11], generalizations of the JKO scheme, in the 
spirit of prediction-correction algorithms (or catching-up algorithms for sweeping 
processes, see for example [12]) might be expected to provide well-posedness results. 

One may also wonder whether it could be possible to recover evolution equations 
for a single species, as in the case of crowd motion models ([11]). In this situation, 
species 2 is replaced by empty space, which can be moved at no cost. Our proposed 
model could be seen as an attempt to replace a unilateral constraint pi < 1 (which 
is quite delicate to handle, see again [11]), by an equality pi + P2 = 1 with p2 > 0. 
It can be done formally by lowering the importance of species 2, more precisely by 
deflning the correction velocity as the one which minimizes a weighted norm: 

w = argmin/ /^^(pi, p2)|v|^, 
Jn 

where the minimum is taken among all those flelds which satisfy V-w = — V-(piUi) 
(in the case U2 = 0), and 

ks{pi,p2) = pi ^£p2' 
One recovers formally an evolution equation for pi with unilateral constraint pi < 1, 
yet with a signiflcant difference: if one considers a saturated zone of 2 surrounded 
by a saturated zone of 1, the asymptotic {e 0) model sees the inclusion as globally 
incompressible, which is not the case if one imposes simply pi < 1. 
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